home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / newmat03.lha / newmat03 / cholesky.cxx < prev    next >
C/C++ Source or Header  |  1993-08-08  |  916b  |  38 lines
  1. //$$ cholesky.cxx                     cholesky decomposition
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.hxx"
  8.  
  9. #include "newmat.hxx"
  10.  
  11.  
  12. /********* Cholesky decomposition of a positive definite matrix *************/
  13.  
  14. // Suppose S is symmetrix and positive definite. Then there exists a unique
  15. // lower triangular matrix L such that L L.t() = S;
  16.  
  17. ReturnMatrix Cholesky(const SymmetricMatrix& S)
  18. {
  19.    int nr = S.Nrows();
  20.    LowerTriangularMatrix T(nr);
  21.    real* s = S.Store(); real* t = T.Store(); real* tk = t;
  22.  
  23.    for (int i=0; i<nr; i++)
  24.    {
  25.       real* tj = t; real* ti = tk;
  26.       for (int j=0; j<=i; j++)
  27.       {
  28.      tk = ti; real sum = 0.0; int k = j;
  29.      while (k--) { sum += *tj++ * *tk++; }
  30.      if (i==j) *tk++ = sqrt(*s++ - sum);
  31.          else *tk++ = (*s++ - sum) / *tj++;
  32.       }
  33.    }
  34.    T.Release();              // not wanted if routine avoids return init
  35.    return T;
  36. }
  37.  
  38.